﻿//--------------------------------------------------------------------------------------------------------------------------------
// <copyright file="SlaeAssembler.cs" company="Exul">
//     Copyright (c) Exul. All rights reserved.
// </copyright>
// <author>Exul</author>
//--------------------------------------------------------------------------------------------------------------------------------

namespace ExulLibrary.FiniteElementMethod.Statements.AxiallySymmetricRectangular2DScalarMaxwell
{
    using System;

    using ExulLibrary.Mathematics.LinearAlgebra;
    using ExulLibrary.Mathematics.LinearAlgebra.Sequential;

    /// <summary>
    ///     Provides method to assembly a system of linear algebraic equations for the task with the specified parameters.
    /// </summary>
    public sealed class SlaeAssembler
    {
        /// <summary>The magnetic permeability of vacuum.</summary>
        public const double VacuumMagneticPermeability = 4 * Math.PI * 1e-7;

        /// <summary>The first stiffness matrix.</summary>
        private static readonly DenseSymmetricMatrix<double> StiffnessMatrixOne =
            new DenseSymmetricMatrix<double>(
                new[]
                {
                    1.0,
                    1.0, 3.0,
                    -1.0, -1.0, 1.0,
                    -1.0, -3.0, 1.0, 3.0
                });

        /// <summary>The second stiffness matrix.</summary>
        private static readonly DenseSymmetricMatrix<double> StiffnessMatrixTwo =
            new DenseSymmetricMatrix<double>(
                new[]
                {
                    2.0,
                    1.0, 2.0,
                    -2.0, -1.0, 2.0,
                    -1.0, -2.0, 1.0, 2.0
                });

        /// <summary>The third stiffness matrix.</summary>
        private static readonly DenseSymmetricMatrix<double> StiffnessMatrixThree =
            new DenseSymmetricMatrix<double>(
                new[]
                {
                    -8.0,
                    0.0, 8.0,
                    -4.0, 0.0, -8.0,
                    0.0, 4.0, 0.0, 8.0
                });

        /// <summary>The first mass matrix.</summary>
        private static readonly DenseSymmetricMatrix<double> MassMatrixOne =
            new DenseSymmetricMatrix<double>(
                new[]
                {
                    2.0,
                    2.0, 6.0,
                    1.0, 1.0, 2.0,
                    1.0, 3.0, 2.0, 6.0
                });

        /// <summary>The second mass matrix.</summary>
        private static readonly DenseSymmetricMatrix<double> MassMatrixTwo =
            new DenseSymmetricMatrix<double>(
                new[]
                {
                    4.0,
                    2.0, 4.0,
                    2.0, 1.0, 4.0,
                    1.0, 2.0, 2.0, 4.0
                });

        /// <summary>The first right-hand vector.</summary>
        private static readonly DenseVector<double> RightOne =
            new DenseVector<double>(new double[] { 1, 2, 1, 2 });

        /// <summary>The second right-hand vector.</summary>
        private static readonly DenseVector<double> RightTwo =
            new DenseVector<double>(new double[] { 1, 1, 1, 1 });

        /// <summary>Gets or sets the rectangular grid.</summary>
        /// <value>The rectangular grid.</value>
        public RectangularGrid2D Grid
        {
            get;
            set;
        }

        /// <summary>Gets or sets the location of material.</summary>
        /// <value>The location of material.</value>
        public RectangularElementLocation2D MaterialLocation
        {
            get;
            set;
        }

        /// <summary>Gets or sets the location of source.</summary>
        /// <value>The source location.</value>
        public RectangularElementLocation2D SourceLocation
        {
            get;
            set;
        }

        /// <summary>Gets or sets the relative magnetic permeability of material.</summary>
        /// <value>The relative magnetic permeability of material.</value>
        public double MaterialRelativeMagneticPermeability
        {
            get;
            set;
        }

        /// <summary>Gets or sets the electric conductivity of material.</summary>
        /// <value>The electric conductivity of material.</value>
        public double MaterialElectricConductivity
        {
            get;
            set;
        }

        /// <summary>Gets or sets the source value.</summary>
        /// <value>The source value.</value>
        public double SourceValue
        {
            get;
            set;
        }

        /// <summary>Gets or sets the time scheme.</summary>
        /// <value>The time scheme.</value>
        public Scheme TimeScheme
        {
            get;
            set;
        }

        /// <summary>Gets or sets the previous solutions.</summary>
        /// <value>The previous solutions.</value>
        public DenseVector<double>[] PreviousSolutions
        {
            get;
            set;
        }

        /// <summary>Assemblies the system.</summary>
        /// <param name="matrix">The matrix of the system.</param>
        /// <param name="rightHandPart">The right-hand part of the system.</param>
        /// <exception cref="ArgumentNullException">
        ///     <paramref name="matrix"/> or <paramref name="rightHandPart"/> is <see langword="null"/>.
        ///     <para>--or--</para>
        ///     <para><see cref="Grid"/> is <see langword="null"/>.</para>
        ///     <para>--or--</para>
        ///     <para><see cref="MaterialLocation"/> is <see langword="null"/>.</para>
        ///     <para>--or--</para>
        ///     <para><see cref="SourceLocation"/> is <see langword="null"/>.</para>
        ///     <para>--or--</para>
        ///     <para><see cref="TimeScheme"/> is <see langword="null"/>.</para>
        ///     <para>--or--</para>
        ///     <para><see cref="PreviousSolutions"/> is <see langword="null"/>.</para>
        /// </exception>
        /// <exception cref="ArgumentException">
        ///     <paramref name="matrix"/> and <paramref name="rightHandPart"/> have incompatible dimensions.
        ///     <para>--or--</para>
        ///     <para><see cref="Grid"/> is an invalid grid.</para>
        ///     <para>--or--</para>
        ///     <para><see cref="PreviousSolutions"/> has length less than <see cref="TimeScheme"/> uses count of layers.</para>
        /// </exception>
        /// <exception cref="ArgumentOutOfRangeException">
        ///     <see cref="MaterialRelativeMagneticPermeability"/> is negative or equal to zero.
        ///     <para>--or--</para>
        ///     <para><see cref="MaterialElectricConductivity"/> is negative or equal to zero.</para>
        ///     <para>--or--</para>
        ///     <para><see cref="MaterialLocation"/> is out of the grid.</para>
        ///     <para>--or--</para>
        ///     <para><see cref="SourceLocation"/> is out of the grid.</para>
        /// </exception>
        public void Assembly(CrcsSparseSymmetricMatrix<double> matrix, DenseVector<double> rightHandPart)
        {
            #if SAFE

            if (matrix == null)
            {
                throw ThrowHelper.ArgumentIsNull("matrix");
            }

            if (rightHandPart == null)
            {
                throw ThrowHelper.ArgumentIsNull("rightHandPart");
            }

            if (matrix.RowsCount != rightHandPart.Length)
            {
                throw ThrowHelper.MatrixVectorDimensionsAreIncompatible(
                    "matrix", matrix.RowsCount, matrix.ColumnsCount, "rightHandPart", rightHandPart.Length);
            }

            if (this.Grid == null)
            {
                throw ThrowHelper.ArgumentIsNull("Grid");
            }

            if (this.MaterialLocation == null)
            {
                throw ThrowHelper.ArgumentIsNull("MaterialLocation");
            }

            if (this.SourceLocation == null)
            {
                throw ThrowHelper.ArgumentIsNull("SourceLocation");
            }

            if (this.TimeScheme == null)
            {
                throw ThrowHelper.ArgumentIsNull("TimeScheme");
            }

            if (this.PreviousSolutions == null)
            {
                throw ThrowHelper.ArgumentIsNull("PreviousSolutions");
            }

            if (this.Grid.AxisOne.Points[0] < 1e-16)
            {
                throw FiniteElementMethod.ThrowHelper.InvalidAxiallySymmetricGrid(
                    "Grid", this.Grid.AxisOne.Points[0], 1e-16);
            }

            string error;
            if (!this.Grid.AxisOne.IsValid(out error))
            {
                throw FiniteElementMethod.ThrowHelper.InvalidGrid("Grid.AxisOne", error);
            }

            if (!this.Grid.AxisTwo.IsValid(out error))
            {
                throw FiniteElementMethod.ThrowHelper.InvalidGrid("Grid.AxisTwo", error);
            }

            if (this.MaterialRelativeMagneticPermeability < 1e-16)
            {
                throw FiniteElementMethod.ThrowHelper.RelativeMagneticPermeabilityIsNegativeOrEqualToZero(
                    "MaterialRelativeMagneticPermeability", this.MaterialRelativeMagneticPermeability);
            }

            if (this.MaterialElectricConductivity < 1e-16)
            {
                throw FiniteElementMethod.ThrowHelper.ElectricConductivityIsNegativeOrEqualToZero(
                    "MaterialElectricConductivity", this.MaterialElectricConductivity);
            }

            if (this.PreviousSolutions.Length < this.TimeScheme.FunctionsCount)
            {
                throw ThrowHelper.CollectionHasNoEnoughElements(
                    "PreviousSolutions", this.PreviousSolutions.Length, this.TimeScheme.FunctionsCount);
            }

            #endif

            // indexes of local element
            var elementIndexes = new int[4];

            // matrix of local element
            var localMatrix = new DenseSymmetricMatrix<double>(4);

            // mass matrix of local element
            var massMatrix = new DenseSymmetricMatrix<double>(4);

            // local right-hand vector from previous solutions (this vector is used to handle time scheme)
            var timeVector = new DenseVector<double>(4);

            // right-hand vector of local element
            var localVector = new DenseVector<double>(4);

            // convert location of the material in space to location in the grid
            var materialLocationInGrid = this.Grid.GetLocationInGrid(this.MaterialLocation);

            // convert location of the source in space to location in the grid
            var sourceLocationInGrid = this.Grid.GetLocationInGrid(this.SourceLocation);

            var axisOne = this.Grid.AxisOne.Points;
            var axisTwo = this.Grid.AxisTwo.Points;

            var axisOnePointsCount = axisOne.Count;
            var axisTwoPointsCount = axisTwo.Count;

            var nodeIndex = 0;
            var topBorder = axisTwo[0];

            for (var rowIndex = 0; rowIndex < axisTwoPointsCount - 1; rowIndex++)
            {
                var bottomBorder = topBorder;
                topBorder = axisTwo[rowIndex + 1];
                var lengthTwo = topBorder - bottomBorder;

                var rightBorder = axisOne[0];

                // a value indicating whether the current location is between bottom and top borders of the material
                var isInsideMaterialByTwo =
                    (materialLocationInGrid.BottomBorderPointIndex <= rowIndex) &&
                    (rowIndex < materialLocationInGrid.TopBorderPointIndex);

                // a value indicating whether the current location is between bottom and top borders of the source
                var isInsideSourceByTwo =
                    (sourceLocationInGrid.BottomBorderPointIndex <= rowIndex) &&
                    (rowIndex < sourceLocationInGrid.TopBorderPointIndex);

                for (var columnIndex = 0; columnIndex < axisOnePointsCount - 1; columnIndex++)
                {
                    var leftBorder = rightBorder;
                    rightBorder = axisOne[columnIndex + 1];
                    var lengthOne = rightBorder - leftBorder;

                    // a value indicating whether the current location is inside the material
                    var isInsideMaterial =
                        isInsideMaterialByTwo &&
                        (materialLocationInGrid.LeftBorderPointIndex <= columnIndex) &&
                        (columnIndex < materialLocationInGrid.RightBorderPointIndex);

                    // a value indicating whether the current location is inside the source
                    var isInsideSource =
                        isInsideSourceByTwo &&
                        (sourceLocationInGrid.LeftBorderPointIndex <= columnIndex) &&
                        (columnIndex < sourceLocationInGrid.RightBorderPointIndex);

                    // calculate indexes of nodes of the current local element
                    elementIndexes[0] = nodeIndex;
                    elementIndexes[1] = nodeIndex + 1;
                    elementIndexes[2] = nodeIndex + axisOnePointsCount;
                    elementIndexes[3] = elementIndexes[2] + 1;

                    // choose magnetic permeability
                    var mu = VacuumMagneticPermeability;
                    if (isInsideMaterial)
                    {
                        mu *= this.MaterialRelativeMagneticPermeability;
                    }

                    // calculate stiffness matrix
                    localMatrix.Assign(StiffnessMatrixOne);

                    var factorOne = lengthOne * lengthOne / (12.0 * lengthTwo * mu);
                    var factorTwo = leftBorder * lengthOne / (6.0 * lengthTwo * mu);
                    localMatrix.Add(factorOne, StiffnessMatrixTwo, factorTwo);

                    localMatrix.Add(StiffnessMatrixThree, lengthTwo / (12.0 * mu));

                    var factorFour = Math.Log(rightBorder / leftBorder) * lengthTwo / (6.0 * lengthOne * lengthOne * mu);
                    localMatrix.Values[0] += factorFour * 2.0 * rightBorder * rightBorder;
                    localMatrix.Values[1] -= factorFour * 2.0 * leftBorder * rightBorder;
                    localMatrix.Values[2] += factorFour * 2.0 * leftBorder * leftBorder;
                    localMatrix.Values[3] += factorFour * rightBorder * rightBorder;
                    localMatrix.Values[4] -= factorFour * leftBorder * rightBorder;
                    localMatrix.Values[5] += factorFour * 2.0 * rightBorder * rightBorder;
                    localMatrix.Values[6] -= factorFour * leftBorder * rightBorder;
                    localMatrix.Values[7] += factorFour * leftBorder * leftBorder;
                    localMatrix.Values[8] -= factorFour * 2.0 * leftBorder * rightBorder;
                    localMatrix.Values[9] += factorFour * 2.0 * leftBorder * leftBorder;

                    if (isInsideMaterial)
                    {
                        timeVector.Clear();

                        // calculate local vector from all previous solutions used in the time scheme
                        var layerIndex = 1;
                        while ((layerIndex < this.TimeScheme.FunctionsCount) && this.TimeScheme.IsFunctionDefined(layerIndex))
                        {
                            var previousSolution = this.PreviousSolutions[layerIndex - 1];
                            var timeDerivative = this.TimeScheme.CalculateDerivativeOnFirstLayer(layerIndex);

                            timeVector.Values[0] -= timeDerivative * previousSolution.Values[elementIndexes[0]];
                            timeVector.Values[1] -= timeDerivative * previousSolution.Values[elementIndexes[1]];
                            timeVector.Values[2] -= timeDerivative * previousSolution.Values[elementIndexes[2]];
                            timeVector.Values[3] -= timeDerivative * previousSolution.Values[elementIndexes[3]];

                            layerIndex++;
                        }

                        // calculate mass matrix
                        massMatrix.Assign(MassMatrixOne);

                        factorOne = lengthOne * lengthOne * lengthTwo * this.MaterialElectricConductivity / 72.0;
                        factorTwo = leftBorder * lengthOne * lengthTwo * this.MaterialElectricConductivity / 36.0;
                        massMatrix.Add(factorOne, MassMatrixTwo, factorTwo);

                        var derivative = this.TimeScheme.CalculateDerivativeOnFirstLayer(0);

                        // add mass matrix to local matrix
                        localMatrix.Add(massMatrix, derivative);

                        // calculate local vector
                        massMatrix.Multiply(timeVector, localVector);

                        rightHandPart.Add(elementIndexes, localVector);
                    }
                    else if (isInsideSource)
                    {
                        localVector.Assign(RightOne);

                        factorOne = this.SourceValue * lengthOne * lengthOne * lengthTwo / 12.0;
                        factorTwo = this.SourceValue * leftBorder * lengthOne * lengthTwo / 4.0;
                        localVector.Add(factorOne, RightTwo, factorTwo);

                        rightHandPart.Add(elementIndexes, localVector);
                    }

                    matrix.Add(elementIndexes, elementIndexes, localMatrix);

                    nodeIndex++;
                }

                nodeIndex++;
            }

            nodeIndex = 0;

            // set up zero conditions on the bottom border
            // elements on the bottom border have only one non-diagonal element, exclude left bottom corner element
            Array.Clear(matrix.TriangleValues, 0, axisOnePointsCount - 1);
            Array.Clear(rightHandPart.Values, 0, axisOnePointsCount);

            while (nodeIndex < axisOnePointsCount)
            {
                matrix.DiagonalValues[nodeIndex] = 1.0;

                nodeIndex++;
            }

            // set up zero conditions on the left and the right borders
            for (var rowIndex = 1; rowIndex < axisTwoPointsCount; rowIndex++)
            {
                // the element on the left border has two non-diagonal elements
                var rowBegin = matrix.Begins[nodeIndex];
                matrix.TriangleValues[rowBegin] = 0.0;
                matrix.TriangleValues[rowBegin + 1] = 0.0;

                matrix.DiagonalValues[nodeIndex] = 1.0;
                rightHandPart.Values[nodeIndex] = 0.0;

                nodeIndex += axisOnePointsCount - 1;

                // the element on the right border has three non-diagonal elements
                rowBegin = matrix.Begins[nodeIndex];
                matrix.TriangleValues[rowBegin] = 0.0;
                matrix.TriangleValues[rowBegin + 1] = 0.0;
                matrix.TriangleValues[rowBegin + 2] = 0.0;

                matrix.DiagonalValues[nodeIndex] = 1.0;
                rightHandPart.Values[nodeIndex] = 0.0;

                nodeIndex++;
            }

            nodeIndex -= axisOnePointsCount;

            // set up zero conditions on the top border
            // elements on the top border have three non-diagonal elements
            Array.Clear(
                matrix.TriangleValues,
                matrix.Begins[nodeIndex],
                3 * (matrix.Begins[nodeIndex + axisOnePointsCount - 1] - matrix.Begins[nodeIndex + 1]));
            Array.Clear(
                rightHandPart.Values,
                nodeIndex + 1,
                axisOnePointsCount - 2);

            nodeIndex++;

            for (var columnIndex = 1; columnIndex < axisOnePointsCount - 1; columnIndex++)
            {
                matrix.DiagonalValues[nodeIndex] = 1.0;

                nodeIndex++;
            }
        }
    }
}